clear all
% load productivities and F distributions (from COUNTERFACTUAL_A_gamma_0.m)

% Simulate (gamma = 1 x estimated value), general equilibrium, revenue
% neutral

global Dividend

Dividend=-688.4+8.3-30;
tau = 0.1;
load fitted_model
% Given Par2 optimal:
n = 50;
r = 0.036;
Par1 = zeros(35,1);
load ParEst_eh0_time0_Mex ParEst_eh0_time0_Mex
load ParEst_eh0_time1_Mex ParEst_eh0_time1_Mex
Par1(1:29) = ParEst_eh0_time0_Mex;
Par1(21) = 10^Par1(21);
Par1(22:27) = Par1(22:27)*1000;
Par1(30) = 500*ParEst_eh0_time1_Mex(1); 
Par1(31) = 500*ParEst_eh0_time1_Mex(2);
Par1(32) = ParEst_eh0_time1_Mex(3);
Par1(33) = ParEst_eh0_time1_Mex(4);
Par1(34) = ParEst_eh0_time1_Mex(5);
Par1(35) = ParEst_eh0_time1_Mex(6);
alphaf = Par1(17); betaf = Par1(18); alphai = Par1(19); betai = Par1(20); 
alphaf_2=Par1(32); betaf_2=Par1(33); alphai_2=Par1(34); betai_2=Par1(35);

load matrix_time1_age0_eh0_comp1.out
% Health transitions and unhealthy proportion
phl_data=matrix_time1_age0_eh0_comp1(1);
plh_data=matrix_time1_age0_eh0_comp1(2);
prob_l=matrix_time1_age0_eh0_comp1(3);

load W_time1_age0_eh0_comp11.out
w_f_g=W_time1_age0_eh0_comp11; % Read wages, g and f distributions
wf_1=w_f_g(:,1);        % Wage in formal sector
wi_1=w_f_g(:,2);        % Wage in informal sector
wf_2=w_f_g(:,7);      % Wage in formal sector
wi_2=w_f_g(:,8);      % Wage in informal sector

wht_f_1=(wf_1(n)-wf_1(1))/(n-1);
wht_i_1=(wi_1(n)-wi_1(1))/(n-1);
wht_f_2=(wf_2(n)-wf_2(1))/(n-1);
wht_i_2=(wi_2(n)-wi_2(1))/(n-1);

p_f_1_est=p_f_1;
p_i_1_est=p_i_1;
p_f_2_est=p_f_2;
p_i_2_est=p_i_2;

% Generate log productivity main pctiles 
[pvector_f_1_est,Fvector_f_1_est]=main_quantiles2(1-Ff_1,p_f_1_est,n);
[pvector_i_1_est,Fvector_i_1_est]=main_quantiles2(1-Fi_1,p_i_1_est,n);
[pvector_f_2_est,Fvector_f_2_est]=main_quantiles2(1-Ff_2,p_f_2_est,n);
[pvector_i_2_est,Fvector_i_2_est]=main_quantiles2(1-Fi_2,p_i_2_est,n);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solve for wage offer parameters

alphaf_old=alphaf;
betaf_old=betaf;
alphai_old=alphai;
betai_old=betai;

iter_wages=1;

for alphaf=alphaf_old*0.99%(0.97:0.01:1.00)%0.99 1 1.01 1 0.99
    for betaf=betaf_old*0.99%(0.98:0.01:1.01)%0.98 0.97 0.97 0.98 0.99
        for alphai=alphai_old*0.99%(0.97:0.01:1.00)%1.01 1 1 1 0.99
            for betai=betai_old*1.05%(1.02:0.01:1.05)%1.02 1.03 1.04 1.05 1.05
                
                Par1(17) = alphaf;
                Par1(18) = betaf;
                Par1(19) = alphai;
                Par1(20) = betai;
                
              [Ff_1,ff_1,Fi_1,fi_1,Ff_2,ff_2,Fi_2,fi_2]=solve_wage_offer_dist(Par1,wf_1,wi_1,wf_2,wi_2,n);
              [Wff_h,Wfi_h,Wfn_h,Wif_h,Wnf_h,Wii_h,Win_h,Wni_h,Wnn_h,Wff_l,Wfi_l,Wfn_l,Wif_l,Wnf_l,Wii_l,Win_l,Wni_l,Wnn_l,Rw_fi_ni_l,Rw_ii_ni_l,Rw_ii_in_l,Rw_if_in_l,...
    Rw_fi_ni_h,Rw_ii_ni_h,Rw_ii_in_h,Rw_if_in_h]=solve_value_function_A1_1(Par1,Ff_1,wf_1,Fi_1,wi_1,Ff_2,wf_2,Fi_2,wi_2,n);
              [HHstates,wage_moments,transitions,phl,plh,gff,gfi,gif,gii,gfn,gnf,gin,gni]=data_simulation_A12_end(Par1,ff_1,wf_1,fi_1,wi_1,ff_2,wf_2,fi_2,wi_2,Wff_h,Wfi_h,Wfn_h,Wif_h,Wnf_h,Wii_h,Win_h,Wni_h,Wnn_h,Wff_l,Wfi_l,Wfn_l,Wif_l,Wnf_l,Wii_l,Win_l,Wni_l,Wnn_l);
                      % marginal distributions
                mff=HHstates(1+9);
                mfi=HHstates(2+9);
                mfn=HHstates(3+9);
                mif=HHstates(4+9);
                mnf=HHstates(5+9);
                mii=HHstates(6+9);
                min_=HHstates(7+9);
                mni=HHstates(8+9);
                mnn=HHstates(9+9);
                mf_1=mff+mfi+mfn;
                mi_1=mif+mii+min_;
                mf_2=mff+mif+mnf;
                mi_2=mfi+mii+mni;
                
                gf_1 = zeros (50,1); gf_2 = zeros (50,1); gi_1 =  zeros (50,1); gi_2 = zeros (50,1);
 for i = 1:50
 gf_1(i) =  (mff*(sum(gff(i,:))) + mfi*(sum(gfi(i,:))) + mfn*gfn(i))/(mff+mfi+mfn);   
 end

  for i = 1:50
 gi_1(i) =  (mif*(sum(gif(i,:))) + mii*(sum(gii(i,:))) + min_*gin(i))/(mif+mii+min_);   
  end
 
   for i = 1:50
 gf_2(i) =  (mff*(sum(gff(:,i))) + mif*(sum(gif(:,i))) + mnf*gnf(i))/(mff+mif+mnf);   
   end
 
    for i = 1:50
 gi_2(i) =  (mfi*(sum(gfi(:,i))) + mii*(sum(gii(:,i))) + mni*gni(i))/(mfi+mii + mni);   
    end
 
                    
                % firm size
                el_f_1=mf_1*gf_1./ff_1;
                % el_f_1(1)=el_f_1(2); % may need smoothing
                % el_f_1(n)=el_f_1(n-1);
                el_i_1=mi_1*gi_1./fi_1;
                % el_i_1(1)=el_i_1(2); % may need smoothing
                % el_i_1(n)=el_i_1(n-1);
                el_f_2=mf_2*gf_2./ff_2;
                % el_f_2(1)=el_f_2(2); % may need smoothing
                % el_f_2(n)=el_f_2(n-1);
                el_i_2=mi_2*gi_2./fi_2;
                % el_i_2(1)=el_i_2(2); % may need smoothing
                % el_i_2(n)=el_i_2(n-1);
                
                der_el_f_1=[Inf;el_f_1(2:n)-el_f_1(1:n-1)]/wht_f_1;
                der_el_i_1=[Inf;el_i_1(2:n)-el_i_1(1:n-1)]/wht_i_1;
                der_el_f_2=[Inf;el_f_2(2:n)-el_f_2(1:n-1)]/wht_f_2;
                der_el_i_2=[Inf;el_i_2(2:n)-el_i_2(1:n-1)]/wht_i_2;
                
                % productivity
                p_f_1=(wf_1+el_f_1./der_el_f_1)*(1+tau);
                p_i_1=wi_1+el_i_1./der_el_i_1;
                p_f_2=(wf_2+el_f_2./der_el_f_2)*(1+tau);
                p_i_2=wi_2+el_i_2./der_el_i_2;
                
                [pvector_f_1,Fvector_f_1]=main_quantiles2(1-Ff_1,p_f_1,n); 
                [pvector_i_1,Fvector_i_1]=main_quantiles2(1-Fi_1,p_i_1,n);
                [pvector_f_2,Fvector_f_2]=main_quantiles2(1-Ff_2,p_f_2,n);
                [pvector_i_2,Fvector_i_2]=main_quantiles2(1-Fi_2,p_i_2,n);
                
                crit_wages(iter_wages)=max([abs(pvector_f_1./pvector_f_1_est-1);abs(pvector_i_1./pvector_i_1_est-1);abs(pvector_f_2./pvector_f_2_est-1);abs(pvector_i_2./pvector_i_2_est-1)]);
                %crit_wages(iter_wages)=sum([(pvector_f_1-pvector_f_1_est).^2;(pvector_i_1-pvector_i_1_est).^2;(pvector_f_2-pvector_f_2_est).^2;(pvector_i_2-pvector_i_2_est).^2]);

                matrix_wages(iter_wages,:)=[alphaf,betaf,alphai,betai];
                iter_wages=iter_wages+1
                
                
            end
        end
    end
end

% Find solution and run for these values
iter_wages
min(crit_wages)
bestproduc=min(find(crit_wages<=min(crit_wages)))

alphaf=matrix_wages(bestproduc,1)
betaf=matrix_wages(bestproduc,2)
alphai=matrix_wages(bestproduc,3)
betai=matrix_wages(bestproduc,4)

Par1(17) = alphaf;
Par1(18) = betaf;
Par1(19) = alphai;
Par1(20) = betai;

 [Ff_1,ff_1,Fi_1,fi_1,Ff_2,ff_2,Fi_2,fi_2]=solve_wage_offer_dist(Par1,wf_1,wi_1,wf_2,wi_2,n);
              [Wff_h,Wfi_h,Wfn_h,Wif_h,Wnf_h,Wii_h,Win_h,Wni_h,Wnn_h,Wff_l,Wfi_l,Wfn_l,Wif_l,Wnf_l,Wii_l,Win_l,Wni_l,Wnn_l,Rw_fi_ni_l,Rw_ii_ni_l,Rw_ii_in_l,Rw_if_in_l,...
    Rw_fi_ni_h,Rw_ii_ni_h,Rw_ii_in_h,Rw_if_in_h]=solve_value_function_A1_1(Par1,Ff_1,wf_1,Fi_1,wi_1,Ff_2,wf_2,Fi_2,wi_2,n);
              [HHstates,wage_moments,transitions,phl,plh,gff,gfi,gif,gii,gfn,gnf,gin,gni]=data_simulation_A12_end(Par1,ff_1,wf_1,fi_1,wi_1,ff_2,wf_2,fi_2,wi_2,Wff_h,Wfi_h,Wfn_h,Wif_h,Wnf_h,Wii_h,Win_h,Wni_h,Wnn_h,Wff_l,Wfi_l,Wfn_l,Wif_l,Wnf_l,Wii_l,Win_l,Wni_l,Wnn_l);
            
wage_moments_gamma_1_GE_rev_neut = wage_moments;
HH_states_gamma_1_GE_rev_neut = HHstates;
save wage_moments_gamma_1_GE_rev_neut wage_moments_gamma_1_GE_rev_neut
save HH_states_gamma_1_GE_rev_neut HH_states_gamma_1_GE_rev_neut


mff=HHstates(1+9);
mfi=HHstates(2+9);
mfn=HHstates(3+9);
mif=HHstates(4+9);
mnf=HHstates(5+9);
mii=HHstates(6+9);
min_=HHstates(7+9);
mni=HHstates(8+9);
mnn=HHstates(9+9);
Un_Head = mnf + mni + mnn; Un_Spouse = mnn + min_ + mfn; Inf_rate = (mnn + min_ + mni + mii);
Agg_stocks_gamma_1_GE_rev_neut= [Inf_rate; Un_Head; Un_Spouse];
save Agg_stocks_gamma_1_GE_rev_neut Agg_stocks_gamma_1_GE_rev_neut
gf_1 = zeros (50,1); gf_2 = zeros (50,1); gi_1 =  zeros (50,1); gi_2 = zeros (50,1);
 for i = 1:50
 gf_1(i) =  (mff*(sum(gff(i,:))) + mfi*(sum(gfi(i,:))) + mfn*gfn(i))/(mff+mfi+mfn);   
 end

  for i = 1:50
 gi_1(i) =  (mif*(sum(gif(i,:))) + mii*(sum(gii(i,:))) + min_*gin(i))/(mif+mii+min_);   
  end
 
   for i = 1:50
 gf_2(i) =  (mff*(sum(gff(:,i))) + mif*(sum(gif(:,i))) + mnf*gnf(i))/(mff+mif+mnf);   
   end
 
    for i = 1:50
 gi_2(i) =  (mfi*(sum(gfi(:,i))) + mii*(sum(gii(:,i))) + mni*gni(i))/(mfi+mii + mni);   
 end

% have to construct the joint probabilities (for ex, gff(wf_1,wf_2)) or integrate over the
% marginal distributions gf_1 and gf_2, whichever is easier..

mff_h=HHstates(1);
mfi_h=HHstates(2);
mfn_h=HHstates(3);
mif_h=HHstates(4);
mnf_h=HHstates(5);
mii_h=HHstates(6);
min_h=HHstates(7);
mni_h=HHstates(8);
mnn_h=HHstates(9);
mff_l=(mff- mff_h*(1-prob_l))/ prob_l;
mfi_l=(mfi- mfi_h*(1-prob_l))/ prob_l;
mfn_l=(mfn- mfn_h*(1-prob_l))/ prob_l;
mif_l=(mif- mif_h*(1-prob_l))/ prob_l;
mnf_l=(mnf- mnf_h*(1-prob_l))/ prob_l;
mii_l=(mii- mii_h*(1-prob_l))/ prob_l;
min_l=(min_- min_h*(1-prob_l))/ prob_l;
mni_l=(mni- mni_h*(1-prob_l))/ prob_l;
mnn_l=(mnn- mnn_h*(1-prob_l))/ prob_l;

Welfare_total_h=r*(sum(sum(mff_h*Wff_h.*gff+mfi_h*Wfi_h.*gfi+mif_h*Wif_h.*gif+mii_h*Wii_h.*gii)))+r*sum(mfn_h*Wfn_h.*gfn+mnf_h*Wnf_h.*gnf+min_h*Win_h.*gin+mni_h*Wni_h.*gni)+r*mnn_h*Wnn_h;
Welfare_total_l=r*(sum(sum(mff_l*Wff_l.*gff+mfi_l*Wfi_l.*gfi+mif_l*Wif_l.*gif+mii_l*Wii_l.*gii)))+r*sum(mfn_l*Wfn_l.*gfn+mnf_l*Wnf_l.*gnf+min_l*Win_l.*gin+mni_l*Wni_l.*gni)+r*mnn_l*Wnn_l;
Welfare_total=prob_l*Welfare_total_l+(1-prob_l)*Welfare_total_h;
Welfare_formal1_h=(r*(sum(sum(mff_h*Wff_h.*gff+mfi_h*Wfi_h.*gfi)))+r*sum(mfn_h*Wfn_h.*gfn))/(mfn_h+mfi_h+mff_h);
Welfare_informal1_h=(r*(sum(sum(mif_h*Wfi_h.*gif+mii_h*Wii_h.*gii)))+r*sum(min_h*Win_h.*gin))/(mif_h+mii_h+min_h);
Welfare_nonemployed1_h=(r*sum(mnf_h*Wnf_h.*gnf+mni_h*Wni_h.*gni)+r*mnn_h*Wnn_h)/(mnf_h+mni_h+mnn_h);
Welfare_formal2_h=(r*(sum(sum(mff_h*Wff_h.*gff+mif_h*Wif_h.*gif)))+r*sum(mnf_h*Wnf_h.*gnf))/(mnf_h+mif_h+mff_h);
Welfare_informal2_h=(r*(sum(sum(mfi_h*Wfi_h.*gfi+mii_h*Wii_h.*gii)))+r*sum(mni_h*Wni_h.*gni))/(mni_h+mii_h+mfi_h);
Welfare_nonemployed2_h=(r*sum(mfn*Wfn_h.*gfn+min_h*Win_h.*gin)+r*mnn_h*Wnn_h)/(min_h+mfn+mnn_h);
Welfare_formal1_l=(r*(sum(sum(mff_l*Wff_l.*gff+mfi_l*Wfi_l.*gfi)))+r*sum(mfn_l*Wfn_l.*gfn))/(mfn_l+mfi_l+mff_l);
Welfare_informal1_l=(r*(sum(sum(mif_l*Wfi_l.*gif+mii_l*Wii_l.*gii)))+r*sum(min_l*Win_l.*gin))/(mif_l+mii_l+min_l);
Welfare_nonemployed1_l=(r*sum(mnf_l*Wnf_l.*gnf+mni_l*Wni_l.*gni)+r*mnn_l*Wnn_l)/(mnf_l+mni_l+mnn_l);
Welfare_formal2_l=(r*(sum(sum(mff_l*Wff_l.*gff+mif_l*Wif_l.*gif)))+r*sum(mnf_l*Wnf_l.*gnf))/(mnf_l+mif_l+mff_l);
Welfare_informal2_l=(r*(sum(sum(mfi_l*Wfi_l.*gfi+mii_l*Wii_l.*gii)))+r*sum(mni_l*Wni_l.*gni))/(mni_l+mii_l+mfi_l);
Welfare_nonemployed2_l=(r*sum(mfn*Wfn_l.*gfn+min_l*Win_l.*gin)+r*mnn_l*Wnn_l)/(min_l+mfn+mnn_l);
Welfare_formal1=prob_l*Welfare_formal1_l+(1-prob_l)*Welfare_formal1_h;
Welfare_informal1=prob_l*Welfare_informal1_l+(1-prob_l)*Welfare_informal1_h;
Welfare_nonemployed1=prob_l*Welfare_nonemployed1_l+(1-prob_l)*Welfare_nonemployed1_h;
Welfare_formal2=prob_l*Welfare_formal2_l+(1-prob_l)*Welfare_formal2_h;
Welfare_informal2=prob_l*Welfare_informal2_l+(1-prob_l)*Welfare_informal2_h;
Welfare_nonemployed2=prob_l*Welfare_nonemployed2_l+(1-prob_l)*Welfare_nonemployed2_h;
Welfare_gamma_1_GE_rev_neut=[Welfare_total;Welfare_total_h;Welfare_total_l;Welfare_formal1;Welfare_informal1;Welfare_nonemployed1;Welfare_formal2;Welfare_informal2;Welfare_nonemployed2]; % IN THOUSANDS, IF NORMALIZED
Budget_balance= -1306.51*Inf_rate-Dividend % the government pays 1306.51 per family/quarter for SP 
save Welfare_gamma_1_GE_rev_neut Welfare_gamma_1_GE_rev_neut

% Generate LOG WAGES by the F distribution 
[wvector_f_1_est,Fvector_f_1_w]=main_quantiles(1-Ff_1,wf_1,n);
[wvector_i_1_est,Fvector_i_1_w]=main_quantiles(1-Fi_1,wi_1,n);
[wvector_f_2_est,Fvector_f_2_w]=main_quantiles(1-Ff_2,wf_2,n);
[wvector_i_2_est,Fvector_i_2_w]=main_quantiles(1-Fi_2,wi_2,n);
mean_wf_1f=log(sum(wf_1.*ff_1));
mean_wi_1f=log(sum(wi_1.*fi_1));
mean_wf_2f=log(sum(wf_2.*ff_2));
mean_wi_2f=log(sum(wi_2.*fi_2));
wage_moments_under_F=[[wvector_f_1_est;mean_wf_1f],[wvector_i_1_est;mean_wi_1f],[wvector_f_2_est;mean_wf_2f],[wvector_i_2_est;mean_wi_2f]];
save wage_moments_under_f_gamma_1_GE_rev_neut wage_moments_under_F

% Transitions 
Trans_head = [transitions(3); transitions(4); transitions(1); transitions(5); transitions(2); transitions(6); transitions(13); transitions(14)];
Trans_spouse = [transitions(9); transitions(10); transitions(7); transitions(11); transitions(8); transitions(12); transitions(15); transitions(16)];
Transitions_gamma_1_GE_rev_neut = [Trans_head,Trans_spouse];
save Transitions_gamma_1_GE_rev_neut Transitions_gamma_1_GE_rev_neut

% Constructing produtivities for simulations

% marginal distributions
mf_1=mff+mfi+mfn;
mi_1=mif+mii+min_;
mf_2=mff+mif+mnf;
mi_2=mfi+mii+mni;


% firm size
el_f_1=mf_1*gf_1./ff_1; 
% el_f_1(1)=el_f_1(2); % may need smoothing
% el_f_1(n)=el_f_1(n-1);
el_i_1=mi_1*gi_1./fi_1;
% el_i_1(1)=el_i_1(2); % may need smoothing
% el_i_1(n)=el_i_1(n-1);
el_f_2=mf_2*gf_2./ff_2; 
% el_f_2(1)=el_f_2(2); % may need smoothing
% el_f_2(n)=el_f_2(n-1);
el_i_2=mi_2*gi_2./fi_2;
% el_i_2(1)=el_i_2(2); % may need smoothing
% el_i_2(n)=el_i_2(n-1);

% M (total employed workforce) in 2005 = 41441076
% Fraction of firms in the formal sector in 2002 includes self-emp. (ENAMIN/INEGI: 23.7%)
% Number of formal firms (IMSS registry, BOSCH AND CAMPOS-VASQUEZ 2005): 800,000
% Number of firms = 800000/0.237
M_N=41441076/(800000/0.237);

pi_f_1=(p_f_1_est-wf_1*(1+tau)).*el_f_1*M_N/0.237;
pi_i_1=(p_i_1_est-wi_1).*el_i_1*M_N/0.763;
pi_f_2=(p_f_2_est-wf_2*(1+tau)).*el_f_2*M_N/0.237;
pi_i_2=(p_i_2_est-wi_2).*el_i_2*M_N/0.763;


profitFormalFirms_1=1./(M_N*mf_1/0.237)*sum(pi_f_1(1:25).*ff_1(1:25)); % per worker 
profitInformalFirms_1=1./(M_N*mi_1/0.763)*sum(pi_i_1(1:18).*fi_1(1:18)); % per worker 
totalprofitFirms_1=mf_1*profitFormalFirms_1+mi_1*profitInformalFirms_1; 
profitFormalFirms_2=1./(M_N*mf_2/0.237)*sum(pi_f_2(1:25).*ff_2(1:25)); % per worker 
profitInformalFirms_2=1./(M_N*mi_2/0.763)*sum(pi_i_2(1:18).*fi_2(1:18)); % per worker 
totalprofitFirms_2=mf_2*profitFormalFirms_2+mi_2*profitInformalFirms_2;  
Profit_gamma_1_GE_rev_neut=[totalprofitFirms_1;profitFormalFirms_1;profitInformalFirms_1;NaN;totalprofitFirms_2;profitFormalFirms_2;profitInformalFirms_2]; % per worker
save Profit_gamma_1_GE_rev_neut Profit_gamma_1_GE_rev_neut

% normalized firm size
[elvector_f_1,Fvector_f_1]=main_quantiles2(1-Ff_1,el_f_1,n)
[elvector_i_1,Fvector_i_1]=main_quantiles2(1-Fi_1,el_i_1,n)
[elvector_f_2,Fvector_f_2]=main_quantiles2(1-Ff_2,el_f_2,n)
[elvector_i_2,Fvector_i_2]=main_quantiles2(1-Fi_2,el_i_2,n)
mean_el_f_1=sum(el_f_1.*ff_1);
mean_el_i_1=sum(el_i_1.*fi_1);
mean_el_f_2=sum(el_f_2.*ff_2);
mean_el_i_2=sum(el_i_2.*fi_2);
Firm_size_gamma_1_GE_rev_neut=[[elvector_f_1,elvector_i_1,elvector_f_2,elvector_i_2];[mean_el_f_1,mean_el_i_1,mean_el_f_2,mean_el_i_2]];
save Firm_size_gamma_1_GE_rev_neut Firm_size_gamma_1_GE_rev_neut


